home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DLLSQ.FOR < prev    next >
Text File  |  1985-11-29  |  6KB  |  242 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DLLSQ
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE LINEAR LEAST SQUARES PROBLEMS, I.E. TO MINIMIZE
  8. C           THE EUCLIDEAN NORM OF B-A*X, WHERE A IS A M BY N MATRIX
  9. C           WITH M NOT LESS THAN N. IN THE SPECIAL CASE M=N SYSTEMS OF
  10. C           LINEAR EQUATIONS MAY BE SOLVED.
  11. C
  12. C        USAGE
  13. C           CALL DLLSQ (A,B,M,N,L,X,IPIV,EPS,IER,AUX)
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           A      - DOUBLE PRECISION M BY N COEFFICIENT MATRIX
  17. C                    (DESTROYED).
  18. C           B      - DOUBLE PRECISION M BY L RIGHT HAND SIDE MATRIX
  19. C                    (DESTROYED).
  20. C           M      - ROW NUMBER OF MATRICES A AND B.
  21. C           N      - COLUMN NUMBER OF MATRIX A, ROW NUMBER OF MATRIX X.
  22. C           L      - COLUMN NUMBER OF MATRICES B AND X.
  23. C           X      - DOUBLE PRECISION N BY L SOLUTION MATRIX.
  24. C           IPIV   - INTEGER OUTPUT VECTOR OF DIMENSION N WHICH
  25. C                    CONTAINS INFORMATIONS ON COLUMN INTERCHANGES
  26. C                    IN MATRIX A. (SEE REMARK NO.3).
  27. C           EPS    - SINGLE PRECISION INPUT PARAMETER WHICH SPECIFIES
  28. C                    A RELATIVE TOLERANCE FOR DETERMINATION OF RANK OF
  29. C                    MATRIX A.
  30. C           IER    - A RESULTING ERROR PARAMETER.
  31. C           AUX    - A DOUBLE PRECISION AUXILIARY STORAGE ARRAY OF
  32. C                    DIMENSION MAX(2*N,L). ON RETURN FIRST L LOCATIONS
  33. C                    OF AUX CONTAIN THE RESULTING LEAST SQUARES.
  34. C
  35. C        REMARKS
  36. C           (1) NO ACTION BESIDES ERROR MESSAGE IER=-2 IN CASE
  37. C               M LESS THAN N.
  38. C           (2) NO ACTION BESIDES ERROR MESSAGE IER=-1 IN CASE
  39. C               OF A ZERO-MATRIX A.
  40. C           (3) IF RANK K OF MATRIX A IS FOUND TO BE LESS THAN N BUT
  41. C               GREATER THAN 0, THE PROCEDURE RETURNS WITH ERROR CODE
  42. C               IER=K INTO CALLING PROGRAM. THE LAST N-K ELEMENTS OF
  43. C               VECTOR IPIV DENOTE THE USELESS COLUMNS IN MATRIX A.
  44. C               THE REMAINING USEFUL COLUMNS FORM A BASE OF MATRIX A.
  45. C           (4) IF THE PROCEDURE WAS SUCCESSFUL, ERROR PARAMETER IER
  46. C               IS SET TO 0.
  47. C
  48. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  49. C           NONE
  50. C
  51. C        METHOD
  52. C           HOUSEHOLDER TRANSFORMATIONS ARE USED TO TRANSFORM MATRIX A
  53. C           TO UPPER TRIANGULAR FORM. AFTER HAVING APPLIED THE SAME
  54. C           TRANSFORMATION TO THE RIGHT HAND SIDE MATRIX B, AN
  55. C           APPROXIMATE SOLUTION OF THE PROBLEM IS COMPUTED BY
  56. C           BACK SUBSTITUTION. FOR REFERENCE, SEE
  57. C           G. GOLUB, NUMERICAL METHODS FOR SOLVING LINEAR LEAST
  58. C           SQUARES PROBLEMS, NUMERISCHE MATHEMATIK, VOL.7,
  59. C           ISS.3 (1965), PP.206-216.
  60. C
  61. C     ..................................................................
  62. C
  63.       SUBROUTINE DLLSQ(A,B,M,N,L,X,IPIV,EPS,IER,AUX)
  64. C
  65.       DIMENSION A(1),B(1),X(1),IPIV(1),AUX(1)
  66.       DOUBLE PRECISION A,B,X,AUX,PIV,H,SIG,BETA,TOL
  67. C
  68. C     ERROR TEST
  69.       IF(M-N)30,1,1
  70. C
  71. C     GENERATION OF INITIAL VECTOR S(K) (K=1,2,...,N) IN STORAGE
  72. C     LOCATIONS AUX(K) (K=1,2,...,N)
  73.     1 PIV=0.D0
  74.       IEND=0
  75.       DO 4 K=1,N
  76.       IPIV(K)=K
  77.       H=0.D0
  78.       IST=IEND+1
  79.       IEND=IEND+M
  80.       DO 2 I=IST,IEND
  81.     2 H=H+A(I)*A(I)
  82.       AUX(K)=H
  83.       IF(H-PIV)4,4,3
  84.     3 PIV=H
  85.       KPIV=K
  86.     4 CONTINUE
  87. C
  88. C     ERROR TEST
  89.       IF(PIV)31,31,5
  90. C
  91. C     DEFINE TOLERANCE FOR CHECKING RANK OF A
  92.     5 SIG=DSQRT(PIV)
  93.       TOL=SIG*ABS(EPS)
  94. C
  95. C
  96. C     DECOMPOSITION LOOP
  97.       LM=L*M
  98.       IST=-M
  99.       DO 21 K=1,N
  100.       IST=IST+M+1
  101.       IEND=IST+M-K
  102.       I=KPIV-K
  103.       IF(I)8,8,6
  104. C
  105. C     INTERCHANGE K-TH COLUMN OF A WITH KPIV-TH IN CASE KPIV.GT.K
  106.     6 H=AUX(K)
  107.       AUX(K)=AUX(KPIV)
  108.       AUX(KPIV)=H
  109.       ID=I*M
  110.       DO 7 I=IST,IEND
  111.       J=I+ID
  112.       H=A(I)
  113.       A(I)=A(J)
  114.     7 A(J)=H
  115. C
  116. C     COMPUTATION OF PARAMETER SIG
  117.     8 IF(K-1)11,11,9
  118.     9 SIG=0.D0
  119.       DO 10 I=IST,IEND
  120.    10 SIG=SIG+A(I)*A(I)
  121.       SIG=DSQRT(SIG)
  122. C
  123. C     TEST ON SINGULARITY
  124.       IF(SIG-TOL)32,32,11
  125. C
  126. C     GENERATE CORRECT SIGN OF PARAMETER SIG
  127.    11 H=A(IST)
  128.       IF(H)12,13,13
  129.    12 SIG=-SIG
  130. C
  131. C     SAVE INTERCHANGE INFORMATION
  132.    13 IPIV(KPIV)=IPIV(K)
  133.       IPIV(K)=KPIV
  134. C
  135. C     GENERATION OF VECTOR UK IN K-TH COLUMN OF MATRIX A AND OF
  136. C     PARAMETER BETA
  137.       BETA=H+SIG
  138.       A(IST)=BETA
  139.       BETA=1.D0/(SIG*BETA)
  140.       J=N+K
  141.       AUX(J)=-SIG
  142.       IF(K-N)14,19,19
  143. C
  144. C     TRANSFORMATION OF MATRIX A
  145.    14 PIV=0.D0
  146.       ID=0
  147.       JST=K+1
  148.       KPIV=JST
  149.       DO 18 J=JST,N
  150.       ID=ID+M
  151.       H=0.D0
  152.       DO 15 I=IST,IEND
  153.       II=I+ID
  154.    15 H=H+A(I)*A(II)
  155.       H=BETA*H
  156.       DO 16 I=IST,IEND
  157.       II=I+ID
  158.    16 A(II)=A(II)-A(I)*H
  159. C
  160. C     UPDATING OF ELEMENT S(J) STORED IN LOCATION AUX(J)
  161.       II=IST+ID
  162.       H=AUX(J)-A(II)*A(II)
  163.       AUX(J)=H
  164.       IF(H-PIV)18,18,17
  165.    17 PIV=H
  166.       KPIV=J
  167.    18 CONTINUE
  168. C
  169. C     TRANSFORMATION OF RIGHT HAND SIDE MATRIX B
  170.    19 DO 21 J=K,LM,M
  171.       H=0.D0
  172.       IEND=J+M-K
  173.       II=IST
  174.       DO 20 I=J,IEND
  175.       H=H+A(II)*B(I)
  176.    20 II=II+1
  177.       H=BETA*H
  178.       II=IST
  179.       DO 21 I=J,IEND
  180.       B(I)=B(I)-A(II)*H
  181.    21 II=II+1
  182. C     END OF DECOMPOSITION LOOP
  183. C
  184. C
  185. C     BACK SUBSTITUTION AND BACK INTERCHANGE
  186.       IER=0
  187.       I=N
  188.       LN=L*N
  189.       PIV=1.D0/AUX(2*N)
  190.       DO 22 K=N,LN,N
  191.       X(K)=PIV*B(I)
  192.    22 I=I+M
  193.       IF(N-1)26,26,23
  194.    23 JST=(N-1)*M+N
  195.       DO 25 J=2,N
  196.       JST=JST-M-1
  197.       K=N+N+1-J
  198.       PIV=1.D0/AUX(K)
  199.       KST=K-N
  200.       ID=IPIV(KST)-KST
  201.       IST=2-J
  202.       DO 25 K=1,L
  203.       H=B(KST)
  204.       IST=IST+N
  205.       IEND=IST+J-2
  206.       II=JST
  207.       DO 24 I=IST,IEND
  208.       II=II+M
  209.    24 H=H-A(II)*X(I)
  210.       I=IST-1
  211.       II=I+ID
  212.       X(I)=X(II)
  213.       X(II)=PIV*H
  214.    25 KST=KST+M
  215. C
  216. C
  217. C     COMPUTATION OF LEAST SQUARES
  218.    26 IST=N+1
  219.       IEND=0
  220.       DO 29 J=1,L
  221.       IEND=IEND+M
  222.       H=0.D0
  223.       IF(M-N)29,29,27
  224.    27 DO 28 I=IST,IEND
  225.    28 H=H+B(I)*B(I)
  226.       IST=IST+M
  227.    29 AUX(J)=H
  228.       RETURN
  229. C
  230. C     ERROR RETURN IN CASE M LESS THAN N
  231.    30 IER=-2
  232.       RETURN
  233. C
  234. C     ERROR RETURN IN CASE OF ZERO-MATRIX A
  235.    31 IER=-1
  236.       RETURN
  237. C
  238. C     ERROR RETURN IN CASE OF RANK OF MATRIX A LESS THAN N
  239.    32 IER=K-1
  240.       RETURN
  241.       END
  242.